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Abstract. Motivated by the recent discovery of a spin liquid phase for the Hubbard 
model on the honeycomb lattice at half- filling |llj . we apply both perturbative and 
non-perturbative techniques to derive effective spin Hamiltonians describing the low- 
energy physics of the Mott-insulating phase of the system. Exact diagonalizations of 
the so-derived models on small clusters are performed, in order to assess the quality 
of the effective low-energy theory in the spin-liquid regime. We show that six-spin 
interactions on the elementary loop of the honeycomb lattice are the dominant sub- 
leading effective couplings. A minimal spin model is shown to reproduce most of the 
energetic properties of the Hubbard model on the honeycomb lattice in its spin-liquid 
phase. Surprisingly, a more elaborate effective low-energy spin model obtained by a 
systematic graph expansion rather disagrees beyond a certain point with the numerical 
results for the Hubbard model at intermediate couplings. 
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1. Introduction 

The quest for exotic quantum phases lacking conventional long-range order in dimensions 
higher than one has recently attracted enormous interest, for fascinating unusual 
properties, such as quantum number fractionalization and/or anyonic excitations, are 
expected to emerge from certain spin Hamiltonians [HE]. In this context, magnetic 
systems are an important playground, since exactly solvable models displaying such 
properties have been recently introduced in the literature pQ. 

Common wisdom has it that, if experimental realizations of such exotic phases are 
to be found, one should look at Hubbard-like systems deep into their strongly interacting 
regime, on frustrated geometries disfavoring more conventional Neel order. Indeed, 
in the limit of large on-site interactions U, where charge fluctuations are inhibited, 
the Hubbard Hamiltonian maps to the Heisenberg model with nearest-neighbor (NN) 
interactions only, for which numerical evidence for a gapped quantum spin liquid (QSL) 
ground-state has been recently found on the highly frustrated kagome lattice [3J. 

An alternative route has been pursued in recent years, by focusing on two- 
dimensional Hubbard models in the regime of intermediate U, where charge fluctuations 
soften the Mott insulating phase. There is indeed convincing evidence that a QSL 
phase is realized on the frustrated triangular lattice [U El [6j U\ El 13 [10]. Sizeable 
charge fluctuations manifest themselves as long-range and/or multi-body effective spin 
interactions, beyond NN spin-exchange. In particular, it has been shown that four- 
spin exchanges are the dominant correction to the Heisenberg model on the triangular 
geometry, accounting for the emergence of a QSL in this case [5j ITU] . 

Whichever of the aforementioned routes is followed, it appears that frustration is an 
essential ingredient for a QSL. In face of this common expectation, the recent discovery 
[TT] of a QSL phase for the Hubbard model on the unfrustrated honeycomb lattice, 
at half-filling and intermediate U, is consequently very surprising and has attracted 
enormous interest. In particular, under the light of our previous discussion, the question 
of what the effective spin interactions are for stabilizing a QSL in this case naturally 
arises. 

An effective spin model comprising terms of up to fourth-order in t/U is basically a 
frustrated J\-Ji model, since four-spin interactions, the dominant corrections to the NN 
Heisenberg model on square and triangular geometries [ID] , are strongly suppressed on 
the honeycomb lattice, for it lacks length-four loops. Due to this fact, it has been argued 
[T2| [13] |H] that the emergence of a QSL can be simply accounted to the frustrating 
next-NN coupling J 2 , a claim that has spurred a number of works on the so-called Jy 
J2 model [15], [161 El [El HDJ- Although such studies differ in detail, they commonly 
detect a quantum critical point, at J%j 3\ ~ 0.2, between long-range ordered Neel and a 
magnetically disordered phase. The precise nature of this non-magnetic state remains 
somehow unclear, albeit most numerical results point to a valence bond solid (VBS) 
[T5J [TBI [13 US]- However, and interestingly enough, in Ref. [JJJ] a variational approach 
finds that such VBS becomes unstable towards a gapped spin liquid if large enough 



Effective Spin Couplings in the Honeycomb Mott Insulator 



3 



length scales are taken into account. 

Independently of the nature of magnetically disordered phase stabilized for the J\- 
J2 model, it is crucial to gauge its validity as a low-energy theory for the honeycomb 
lattice Hubbard Hamiltonian in the intermediate range of t/U, where the QSL emerges. 
Indeed, evidence indicating a more involved situation has been found in Ref. [2U]: the 
bare series expansion (of up to order 14 th ) is not converged for values of t/U leading 
to the QSL. Furthermore, a non-perturbative derivation of effective spin interactions, 
obtained via graph-based continuous unitary transformations [20] , yields a relatively 
small ratio J2/J1 < 0.2, insufficient to destabilize the long-range ordered Neel phase. It 
is thus clear that a more thorough analysis is called for, this being our main goal in the 
present work. 

We employ three state-of-the-art methods to derive effective spin models for the 
Hubbard model on the honeycomb lattice at half-filling, in order to identify dominant 
effective interactions as a function of t/U. Interestingly, among the vast number of 
effective spin couplings shooting up in the vicinity of the quantum phase transition to 
the semi-metal, we find that the largest correction to the NN Heisenberg interaction are 
six-spin interactions located on single hexagons; the frustrating next-NN exchange J2 is 
considerably smaller. We assess the accuracy of the effective low-energy theory, for t/U 
yielding the QSL, by performing exact diagonalizations (EDs) on small clusters. 

We organize the paper as follows. In Sec.|2]we write down the Hubbard Hamiltonian 
and a generic effective spin model, accordingly fixing the notation employed throughout 
the paper. The methods used in deriving the low-energy spin model, namely perturbative 
and graph-based continuous unitary transformations (respectively, pCUTs and gCUTs), 
as well as the contractor renormalization (CORE) group, are described in Sec. [3j After 
comparing the outcomes from both approaches in Sec. |3.3[ we present results from EDs 
in Sec. |4| Finally, we summarize our results in Sec. [5] 



We consider the single-band Hubbard model studied in Ref. [UJ , defined on a honeycomb 
geometry and at half-filling, that reads 



= c ia Ci a is the occupation number operator for fermions with spin a at the site 
i of the honeycomb lattice, and t is the amplitude for hoppings taking place between 
NN sites, on this lattice. Quantum Monte Carlo (QMC) simulations pj] show 

that a QSL phase is stabilized for moderately strong couplings, 0.233 < t/U < 0.286 
(4.3 > U/t > 3.9), and it is our main purpose here to compute effective spin interactions 
for Eq. Q in this range. 

Due to the SU(2) symmetry of the Hubbard Hamiltonian, a generic effective model 
for its Mott insulating phase can be expressed in terms of products of spin-half operators 



2. Model 




(1) 
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-f^spin = E + ^ Jij (S{ ■ Sj^j + ^ Kij k i ySi ■ Sj^j (^Sk ■ Si^j (2) 

+ Lijki mn (Si • Sj^j (j^k ■ Si^j (j~> m • S^j + . . . , 

i,j,k,l,m,n 

where E denotes a constant energy shift. J^, K^ki and Lijkimn respectively denote 
coupling constants for the various two-, four-, and six-spin exchanges and the . . . refer 
to analogue expressions involving more than six spins. Unfortunately, however, the 
situation is complicated by the fact that multi-spin exchanges involving eight or more 
spins form an over-complete set of operators [2T]. As a consequence, no unique solution 
in terms of spin operators can be obtained. While the effective Hamiltonian remains 
well defined in terms of matrix elements in the spin basis, it is only the algebraic 
representation of the effective Hamiltonian in terms of spin operators that is not unique. 

In the next section we will describe different approaches to derive such effective 
low-energy models. Furthermore, we extract and discuss the most important corrections 
to the nearest-neighbor Heisenberg exchange giving rise to a minimal magnetic model 
capturing the essential physics of the full Mott phase. This minimal model as well as a 
more elaborate effective spin model is then analysed afterwards. 

3. Methods and Effective Models 

In this section, we discuss details concerning the numerical methods employed in 
obtaining the effective spin interactions for the Hubbard model, Eq. 0, at half- 
filling and in the regime of strong to intermediate couplings. We start by the so- 
called continuous unitary transformations which allow us to gain a global view on the 
effective spin model and its most important spin interactions. Afterwards, we apply 
the contractor renormalization technique to confirm the behaviour of the dominant 
effective spin couplings. We find that both methods essentially agree when considering 
the same minimal set of considered clusters in both calculations. The latter motivates 
the definition of a minimal magnetic model. 

3.1. CUTs and full effective spin model 

We have calculated the dependence of the magnetic exchange couplings on t/U using 
perturbative continuous unitary transformations (pCUTs) [22l |23l El] along the lines of 
Ref. pm |25]. The pCUT provides the magnetic couplings as series expansions in t/U. 
Note that the spectrum of the Hubbard model is symmetric at half filling under the 
exchange t •f-)- —t and therefore only even order contributions are present [26] . We have 
determined all two-spin, four-spin and six-spin interactions up to order 14. In order 14 
one has to calculate 345 topologically different graphs. 

The bare series do not converge in the spin-liquid region of the Mott phase [20] . 
This is different to the recently analysed case of the Hubbard model on the triangular 
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lattice pm [20] where the intermediate non-magnetic phase is already realized for smaller 
ratios of the bandwidth W to the interaction U. It is therefore mandatory to apply 
resummation schemes of the series which is very complicate to be performed reliably for 
the full set of exchange couplings. We have therefore applied the recently developped 
graph-based continuous unitary transformations (gCUTs) [20] in order to derive non- 
perturbatively the effective spin model. In the following we use the gCUT results to 
study the physics of the full low-energy spin model. The pCUT results are only used to 
analyze the most important terms in the effective model. 

The general properties of the gCUT are discussed in Ref . [20] . Here we only mention 
the specific approach for the current problem. The basic idea is to use a CUT to map 
H unitarily to an effective Hamiltonian H c g which has the special property that the 
block without double occupancies representing the effective low-energy spin model is 
decoupled from the rest of the Hamiltonian. 

In the gCUT one generates all topologically distinct connected graphs Q v of the 
lattice and one sorts them by their number of sites n. On each graph Q v a CUT 
is implemented by setting up the finite number of flow equations [27J [28j [29] . A 
continuous auxiliary variable I is introduced defining the /-dependent Hamiltonian 
H G "(l):=U\l)H Gv U(l). Then the flow equation is given by 



where rj(l):=—W (l)(diU(l)) is the anti-Hermitian generator of the unitary transforma- 
tion. At the end of the flow I = oo one obtains an effective graph-dependent Hamiltonian 



number of double occupancies and H denote all remaining non-diagonal processes. We 
stop the flow for each graph once the residual off-diagonality (ROD) is below 10~ 9 . Here 
the ROD is defined as the square root of the sum of all non-diagonal elements which 
connect to the low-energy subspace. Finally, one obtains the effective spin model H sp i n in 
the thermodynamic limit by substraction of subcluster contributions and by embedding 
back the graphs on the honeycomb lattice [20j . 

The numerical effort scales exponentially with the number of sites n of a given 
graph. Here we have treated all graphs up to 7 sites. This amounts to solving up to 
one million differential equations for the most demanding symmetry sector of a single 
graph for each value of t/U. Physically, one captures all charge fluctuations on this 
length scale. Let us remark that it is of course possible to only include a restricted set 
of graphs in the gCUT calculation. This will be used below for a direct comparison with 
CORE. Finally, one obtains for a given t/U the effective spin model non-perturbatively. 

A view on the full effective spin model obtained by gCUTs for different values of 
t/U is shown in Fig. IT] Let us remark that for gCUT(7) the full low-energy model can be 



^(0 = ^(0,^(01 



(3) 





Figure 1. (Color online) Coefficients of all gCUT(7) terms in the effective Hamiltonian 
for different U/t. Terms with different number of spin operators 0, 2, 4, or 6 are 
separated by vertical dashed lines. Each block containing the same number of spin 
operators is sorted by the size of the coefficient. 



uniquely expressed in terms of spin operators, because multi-spin interactions involving 
more than six sites are not yet allowed for this cluster size. 

There are several implications which can be directly read off from Fig. [T} First, 
only at rather small values of t/U < 0.1 a clear hierarchy in the amplitudes of the 
effective spin operators can be seen. This corresponds to the purely perturbative regime 
where the importance of terms is fully given by the perturbative order where a spin 
coupling appears for the first time. In this regime the frustrated next-nearest neighbor 
Heisenberg exchange J2 (order 4 perturbation theory) is the leading correction to the 
dominant Heisenberg interaction J\ (order 2 perturbation theory). 

Second, the situation changes drastically in the intermediate coupling regime where 
the spin-liquid phase is realized. Here, no clear hierarchy in the amplitudes of the 
effective spin couplings can be detected. One is rather confronted with a proliferating 
number of terms which is likely a consequence of the fact that the metal-insulator 
transition is second order. Consequently, the charge gap closes continuously and one 
expects an increasing number of effective spin couplings on increasing length scales 
when approaching the quantum critical point. Let us note that this is different for the 
Hubbard model on the triangular lattice pU]. In this case the metal-insulator transition 
is first order and one has no diverging length scale upon approaching the Mott transition 
from the insulating side. 

Third, and despite the complexity of the full effective spin model, one still expects 
that most of the properties of the Mott phase should be contained in a minimal 
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Figure 2. (Color online) Dominant effective spin interactions for the Hubbard model 
on the honeycomb lattice at half-filling [Eq. Q]: NN and next-NN two-spin couplings 
(Ji and J2, respectively), and six-spin exchange terms with couplings L%, L2, • • • , L5. 

magnetic model where one only considers the most important spin couplings having 
the largest amplitudes. Some evidences for this strategy are given in the next section by 
analysing such a minimal magnetic model. Interestingly, we find that the perturbative 
hierarchy is not valid anymore for intermediate couplings: the most important correction 
to Ji is not J2 but rather six-spin interactions located on hexagons. These six-spin 
interactions arise in order six perturbation theory and they originate dominatly from 
ring exchange processes around the shortest loop of even length on the honeycomb lattice 
corresponding to hexagons. Let us remark that this is very similar to the Hubbard model 
on the triangular lattice |10| . Here the most elementary loop of even length is a four- 
site plaquette. Consequently, four-spin interactions on the elementary plaquette are the 
dominant correction to the Heisenberg model. 

In the remainder of this paper we are aiming at a minimal magnetic model 
which contains only the most important spin interactions in the regime of strong to 
intermediate coupling. In this parameter regime, the minimal model should display the 
low-energy properties of the Hubbard model on the honeycomb lattice. Additionally, we 
expect that going beyond the minimal model by considering the full graph decomposition 
used by gCUT (or alternatively by CORE) yields successfully even better results. The 
most important effective exchanges are depicted in Fig. g NN ( Jf) and next-NN ( J 2 ) 
two-spin interactions appearing in the J1-J2 model, and the six-spin terms with couplings 
Li, L 2 , ■ ■ ■ , L5. The corresponding gCUT and pCUT amplitudes for these couplings are 
displayed in Fig. [3j 

One clearly sees that the bare series of the six-spin interactions are not converged 
in the interesting regime of intermediate t/U confirming the conclusion obtained in 
Ref. J2U]. In contrast, the results obtained by gCUT and self-similar extrapolation of 
the pCUT series are in fair agreement for all displayed couplings. 

In the following we want to further strengthen these findings by using CORE as an 
alternative tool to derive effective low-energy models. 




Figure 3. (Color online) The relative gCUT and pCUT amplitudes J2/J1 (left) and 
Li/Ji (right) as a function of t/U. Thin solid lines refer to the bare pCUT series. Black 
solid lines correspond to self-similar extrapolants as discussed in Rcf. [20] . Symbols 
denote results obtained by gCUTs. 



3.2. Contractor renormalization and minimal model 

It is also possible to derive non-perturbative effective theories by relying on the 
CORE technique [30]. CORE is a real-space renormalization procedure, applicable 
to generic lattice models, where effective models are derived by truncating local degrees- 
of-freedom. Such construction requires the exact computation of low-lying eigenstates 
on finite connected graphs, similarly to what happens in the previously discussed gCUT 
procedure; implementation details are discussed in the literature [301 [31]. The resulting 
effective Hamiltonian is given as a cluster expansion 

^eff = K > ( 5 ) 

9 

where the sum takes place over a set of graphs g, and h c g corresponds to the connected 
term that is obtained by substracting contributions from embedded sub-clusters [3U1 |3T] • 
Such expansion must naturally be truncated at a certain maximum range [32J, but its 
accuracy can in principle be controlled by analyzing the convergence of terms in the 
effective model with increasing range. 

Since it is our goal here to obtain effective spin interactions for the Hubbard model 
[Eq. p])] at half-filling, we select single sites as elementary blocks in the CORE expansion 
[30] l3l] and retain only singly-occupied states on all sites. 

The clusters used in the CORE calculation are shown in Fig. |4} Here two setups 
are considered: i) The first choice of clusters does not contain the four-site star graph 
displayed in Fig. Hie). It therefore corresponds to the leading term of a full graph 
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Figure 4. Clusters considered in the CORE expansion, ordered according to the 
maximum range for effective interactions |32j and comprising: (a) two, (b) three, (c) 
four and (d) six sites. 

expansion in terms of hexagons, ii) The second choice contains additionally the four- 
site star cluster [Fig. 121(c)]. In the hexagon expansion mentioned before such a graph 
would be included only on the level of three-hexagon clusters. Later we see that gCUT 
(restricted to exactly the same graphs) yields almost identical results for both choices 
of clusters. The second choice can then be regarded as an intermediate step between 
the minimal one-hexagon calculation and the full gCUT(7) calculation which includes 
all graphs shown in Fig. E] CIS db subset. 

We start by discussing the first choice. The lowest-order contribution to the NN 
spin-exchange J\ is simply computed by solving the Hubbard model on a two-site cluster 
[Fig. |4](a)]: in this case only eigen-energies are required, since J\ exactly corresponds to 
the singlet-triplet gap: 

Although trivially obtained, this result is non-perturbative, reducing to the well-known 
limit Jx ~ At 2 /U only when t/U < 1 [see Fig. plb)]. Let us mention that the same 
non-perturbative result is obtained when calculating J{ with gCUT(2). 

Longer-range couplings, along with corrections to shorter-ranged ones, are 
obtainable from the analysis of the larger clusters depicted in Fig. |4](b-d). In particular, 
in Fig. [5|a) we plot bare effective couplings, obtained by considering the six-site cluster 
in Fig. |4](d), as a function of t/U . Here the notion 'bare couplings' corresponds to the 
effective exchange couplings of a single cluster without considering substractions and 
embeddings. We first focus on the next-NN coupling J 2 and notice that, particularly in 
the range 0.233 < t/U < 0.286 ((4.3 > U/t > 3.9)) where a QSL appears [llj, J 2 < 0.1 Ji 
and therefore this frustrating interaction is not large enough [151 HS1 El CHI EES] to 
account for the absence of long-range Neel order before a semi-metal phase is stabilized 
|20j . However, we also notice that, remarkably, much larger magnitudes are associated 
to the six-spin terms with couplings Li,L 2 , • • • ,L$ depicted in Fig. [2] (we also briefly 
remark that the signs for their amplitudes are in agreement with those obtained from a 
decomposition of the cyclic 6-spin permutation, cf. Ref. [33]). Furthermore, we observe 
that the amplitudes for longer-range two-spin, as well as four-spin, terms are found to be 
negligible in comparison to J\ and L±, L 2 , ■ ■ ■ , L 5 [34J. One is thus inclined to conclude 




Figure 5. (Color online) (a) Bare amplitudes for of all terms in the effective 
Hamiltonian for Eq. ([IJ as a function of t/U, as obtained from the CORE procedure 
applied to the six-site cluster depicted in Fig. |4]jd). For clarity sake, only dominant 
couplings are shown, following the notation from Fig. [2j all remaining terms have very 
small amplitudes and lie in the filled region at the bottom-left corner, (b-c) Results 
for the two-spin exchanges J\ and J2, as a function of t/U, obtained from indicated 
clusters (we label the results according to the number of sites in a given cluster; see 
Fig. EJ). For each cluster, contributions from embedded clusters comprising lesser sites 
are subtracted, following the standard CORE recipe [30l 131] , In both panels, the 
shaded region corresponds to the QSL phase [TT] . 



that such six-spin interactions are the main ingredients in stabilizing a QSL in the phase 
diagram for Eq. ([T|, a conjecture we aim at testing in Sec. |ij 

Let us mention that such six-spin interactions on local hexagons are generated in 
order six perturbation theory. In this order the couplings L\-L^ have the same absolute 
prefactor with an alternating sign which exactly corresponds to the six-spin interactions 
contained in the six-spin ring exchange permutation operator. But we stress that the 
latter operator also contains two-spin and four-spin interactions of the same magnitude 
which one has to contrast with our finding, both perturbatively and non-perturbatively, 
of suppressed two- and four-spin interactions. The full non-perturbative contribution 
of a single hexagon is therefore dominated by the six-interactions having the largest 
exchange couplings. Finally, we note that the nonperturbatively obtained couplings 
L1-L5 have different prefactors as can be seen in Fig. |5j 

We proceed by analyzing the dependence of the effective couplings on the maximum 
range of the clusters. In order to do so, for each cluster depicted in Fig. [| we subtract 
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connected contributions from shorter-ranged, embedded, sub-clusters, according to the 
standard CORE recipe [301 EE]- We plot results for the two-spin exchanges with 
couplings J\ and J 2 as a function of t/U respectively in Figs. [5](b) and (c), and find 
a seemingly fast convergence with increasing cluster range up to the six-site cluster 
depicted in Fig. Eld). Furthermore, the data in Figs. [5^b) and (c) is in good agreement 



with gCUT results, as we discuss below in Sec. |3.3 

Finally, we comment on the fact that a few clusters larger than the ones in Fig. |4| 
such as the one comprising ten sites and forming an edge-sharing double hexagon, are 
amenable to numerical diagonalization, in principle implying that our CORE expansion 
could be extended to even longer ranges than considered here. To this end a full graph 
decomposition must be implemented for the CORE approach which will be left as a 
task for future research. Here we would rather like to focus on the most minimal setup 
to describe the low-energy physics of the Hubbard model. 

3.3. Minimal effective model 



We compare now the effective models obtained from gCUT (Sec. 3.1) and CORE 



(Sec. 3.2). We first remark that both approaches yield qualitatively similar effective 
Hamiltonians, that display as dominant terms the two- and six-spin exchanges depicted 
in Fig. [2] (see Figs. [3] and [5]). We therefore focus on such terms, that may be regarded 
as defining a minimal effective model, and plot their amplitudes, as a function of t/U in 
Fig. [6j Let us mention that we do not show the constant term of the effective models 
which also displays a t/U dependence and which one therefore has to keep in the exact 
diagonalization when comparing to the Hubbard model. Here three different sets of 
graphs are considered: i) the leading graphs of a hexagon expansion not including the 
four-site star graph (CORE and gCUT), ii) the leading graphs of a hexagon expansion 
including the four-site star graph (CORE and gCUT), and iii) all graphs comprising up 
to seven sites (gCUT). 

We notice that results from both methods are in good quantitative agreement up 
to couplings stabilizing the QSL phase [UJ and beyond the perturbative regime [20] 
for all three choices. Especially, CORE and gCUT yield almost identical results for 
the restricted graph sets i) and ii). This is likely a consequence of the fact that the 
structure of the low-energy subspace is rather constraint due to the high symmetry of 
the considered clusters (the four-site graph as well as the six-site graph has a rotational 
symmetry) plus the number of different spin couplings defined on these graphs is rather 
small. 

The most noticeable effect when going from the minimal choice to the full graph 
decomposition including all graphs up to seven sites is an increasing value for the next- 
NN two-spin exchange J 2- Nevertheless, we also remark all three sets agree in that, 
unlike what is assumed in Refs. [T2], [T3], the effective values for the next-NN two-spin 
exchange J 2 are not large enough to destroy Neel order [T5J [TBI E0 EEH1 EE] , suggesting 
that the QSL phase may instead be accounted for by the six-spin ring exchanges with 
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□ — □ CORE 3 graphs 
gCUT 3 graphs 
□— □ CORE 4 graphs 
gCUT 4 graphs 
C gCUT(7) 



Figure 6. (Color online) Comparison between effective couplings, obtained from 



gCUT (Sec. 3.1) and CORE (Sec. 3.2), for the dominant spin-exchanges depicted in 
Fig. [2} The range of values for t/U leading to the QSL phase, according to Ref. [IT] , 
is highlighted. 



amplitudes Li, L2, ■ ■ ■ L5 in Fig. [2] Furthermore, we show in the next section that the 
minimal effective spin model having the lowest J2 give the best agreement with the 
QMC results at intermediate couplings. 



4. Exact diagonalizations 

We turn now our attention to the characterization of the effective spin model derived 
in the previous sections, either with gCUT method or CORE algorithm. Our goal is 
to show that a reasonably simple and compact effective spin model can describe rather 
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well the energetics of the Hubbard model. In order to do so, we will make a systematic 
comparison of the low-energy properties, and of the local quantities such as double 
occupancy or bond kinetic energy, that we compute using the exact diagonalization 
(ED) technique. More precisely, we have used a standard Lanczos algorithm to obtain 
the ground-state energy for various clusters at fixed total S z , and we have made use of 
all space symmetries (translations and point-group, when available). 



4-1- Study of the minimal effective model 

We start by discussing the properties of the minimal effective spin model obtained for 
the minimal set of graphs not including the four-site star graph. We stress again that 
CORE and gCUT give essentially identical low-energy spin models for this choice of 
graphs. In the following we use the effective model obtained via the CORE algorithm 



(Sec. 3.2) 



We start by comparing the low-energy spectra for the Hubbard [Eq. ([T|] and the 
effective (Sec. [3]) models onaniV = 18 cluster for two values of t/U : t/U = 0.05 (deep 
into the Neel phase) and t/U = 0.25 (QSL phase, according to Ref. In Fig. uta-b), 

we plot so-called Anderson tower of states for the Hubbard model. As we only make use 
of S z symmetry (instead of the total spin), we plot the energy states versus S Z (S Z + 1), 
so that degenerate states with identical energies correspond to spin multiplets. We have 
also computed the one-particle gap, A 1P = (E (N/2 + 1) + E (N/2 - 1) - 2E (N/2))/2, 
where Eo(Nf) is the ground-state energy with N f fermions. For t/U = 0.05 [Fig. ga)], 
the so-called Anderson tower of states is typical of a collinear Neel state [35J with a 
lowest excitation increasing as S(S + 1) (with a slope that should decrease as l/N), 
followed by spin- wave excitations (with energy scale 1/y/N). Moreover, on the scale 
of the figure, all states correspond to spin excitations since A 1P is quite large [36J. 
The same behaviour is reproduced qualitatively and quantitatively by the effective spin 
model, as expected and shown in Fig. [7|c). Note that in the case of the effective spin 
model, we have been able to label each eigenstate with its total spin value S so that we 
plot levels function of S(S + 1). 

When increasing t/U to 0.25, a distinct picture emerges [Fig. [Tj^b)] . While the first 
excitations still correspond to S — 1 and S = 2 states at the T point, there is no longer a 
straight S(S + 1) line of excitations and, even more importantly, the gap A 1P to charge 
excitations is substantially lowered, beyond which a description solely based on spin 
degrees-of-freedom is no longer valid. However, according to QMC results [11] , for this 
coupling the one-particle gap should remain finite in the thermodynamic limit, so that 
our effective spin model could still describe the physics at the lowest energies: indeed, 
its low-energy spectrum, plotted in Fig. fftd), reproduces the lowest spin excitations 
quite accurately. Furthermore, while the quantum numbers seen in the tower-of-states 
are consistent with standard collinear magnetic ordering, the lowest spin excitations do 
not display a clear S(S + 1) behaviour, which might be an indication that there is no 
magnetic ordering. 
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Figure 7. (Color online) (a-b) Low-energy spectrum of the Hubbard model vs 
S Z (S Z + 1) obtained for t/U = 0.05 and t/U = 0.25 on N = 18 honeycomb cluster. 
Different symbols correspond to different points in the Brillouin zone: T point at the 
center, 6-fold degenerate A point and 2-fold degenerate K point at the corners (see 
inset). We also give twice the one-particle gap Aip in each case, (c-d) Low-energy 
spectrum of the CORE effective model vs S(S + 1) obtained on the same N = 18 
honeycomb cluster. 



It would also be desirable to compute correlators for the ground-state of the effective 
model, but this is quite involved since one would also need to renormalize the operators 
when deriving an effective Hamiltonian [301 EH HOI [20] • Nevertheless, some observables 
may straightforwardly be obtained from the ground-state energy per site, eo, by relying 
on the Feynman-Hellmann theorem. For instance, the electronic double occupancy per 
site can be computed from ED of the effective model involving only spin degrees-of- 
freedom, as it has been done in Ref. [TO] using : 

<W> = ^ (7) 

Results are shown in Fig. [8(a) and compared both to the exact result of the Hubbard 
model on N = 18 cluster and to QMC data [37]. Since this is a local quantity, its finite 
size effects are rather small as expected. We observe a good agreement not only for small 
t/U (deep in the Neel phase), but also in the interesting regime t/U ~ 0.25 where QSL 
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is expected to emerge. As a side remark, let us emphasize that the discrepancy between 
results on the Hubbard model and with the CORE effective one on N = 18 can be 
attributed to the role of short-length loops (of length 6) that exist on this small cluster, 
and are not captured by the effective model (see Ref. [10J [15] for a similar discussion). 
In fact, effective model results agree quite well with the data obtained on a much larger 
lattice with QMC. 

In Ref. [TT], the analysis of the behaviour of the kinetic energy density E^ in = 
(—t a( c la c jv + h.c.))/N versus U was also discussed in the context of the formation 
of local moments at strong U/t, in contrast with the itinerant regime for small U/t, 
whereas the QSL is stabilized in between. Using Feynman-Hellmann theorem again, we 
can simply compute this quantity as 

We have chosen a grid of U/t going from 3 to 8 by step of 0.1 and computed the ground- 
state energy for various clusters. Using a finite-difference approximation, we obtain the 
derivative of the kinetic energy that we plot against t/U in Fig. |8](b). Our first comment 
is that we have some agreement with the exact QMC data [Hj that show a maximum 
around 0.14 for U/t ~ 5.0, although this quantity is expected to be quite sensitive to 
details since it is a second derivative of the ground-state energy. Note that QMC data 
are also quite noisy in this large t/U region, and we refer to the original publication [11] 
for additional data : the point is that a maximum of order 0.15 is reached for t/U ~ 0.2. 
Our second remark is that we observe an anomaly around t/U = 0.3 which could signal 
that our CORE truncation is not safe beyond this value. This gives us confidence that 
at least for t/U ~ 0.25, our approach could make sense, and this includes the spin liquid 
phase. Last, since we observe this change of curvature around U/t = 4.5 in a spin-only 
model, it seems to us that the virtual charge fluctuations giving rise to our effective 
low-energy model are sufficient to account for this phenomenon. Additionally, it might 
be not surprising that we underestimate this effect by our minimal magnetic model, 
because we expect that the large number of neglected additional spin operators as well 
as small renormalizations of the treated spin couplings could well lead to an enhanced 
itineracy of the ground state at intermediate t/U values. 

As a conclusion on this part, while we are clearly not able to reach system lengths 
where QSL behavior is expected to occur - which would require several hundred sites - 
our simple effective model is able to reproduce both the low-energy properties and the 
ground-state local properties. It suggests that these properties may not be linked to the 
presence of a QSL and are robust short-distance features. 

4-2. Study beyond the minimal model 

The results of the last subsections are very promising. CORE and gCUT give an almost 
identical minimal effective spin model which is obtained from the single hexagon graph 
and which compares well to the results obtained by QMC for the Hubbard model on 
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Figure 8. (Color online) (a) Double occupancy per site obtained obtained with the 
minimal CORE effective model (3 graphs) on finite clusters with N = 18 to N = 30. 
(b) Derivative of the kinetic energy, as a function of t/U, obtained from the ground- 
state energy of the effective model [Eq. (J8|] on the same clusters. In both cases, 
comparison is made to ED (N = 18) and QMC (N = 72) [37] data of the Hubbard 
model. 

the honeycomb lattice. 

In the following we want to go beyond this minimal choice of graphs and we would 
like to see whether the above findings are confirmed and strengthened. We therefore 
compare the above results for the minimal model to exact diagonalizations for the 
effective spin models obtained i) by including the four-site star graph to the minimal 
set and ii) by including all graphs up to seven sites. For case i) we again use results 
from the CORE algorithm, but we stress that the gCUT on the restricted set of graphs 
gives basically the same low-energy model. For the second case the results of gCUT(7) 
are used. 

We start by comparing the ground-state energy per site as a function of t/U which 
is displayed in Fig. [9] It can be clearly seen that all three different sets of graphs give 
a very similar energy per site in the full Mott phase which also agrees well with results 
for the Hubbard model. This is promising, but it is also a consequence of the fact that 
the ground-state energy per site is a rather unsensitive quantity. 

We therefore turn to more sensitive quantities, namely the double occupancy per 
site and the derivative of the kinetic energy as discussed in the last subsection for 
the minimal effective model. Consequently, details of the considered effective spin 
models matter which is most clearly seen for the derivative of the kinetic energy. 
Interestingly and surprisingly, the agreement between results from the effective spin 
model and results for the Hubbard model gets worst at intermediate couplings. Indeed, 
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Figure 9. (Color online) Ground-state energy per site as a function of t/U. A 
comparison is made to ED (N = 18) and QMC (N = 72) [37] data of the Hubbard 
model. 



already the minimal set of graphs plus the four-site star graph yields the wrong tendency. 
Furthermore, exact diagonalization of the full effective spin model from gCUT(7) does 
not even display the maximum in the derivative of the kinetic energy. 

Altogether, we find that the inclusion of more graphs in our calculation gives a 
reduced agreement at intermediate couplings in contrast to our expectation. This trend 
is most probably neither a numerical problem nor a difference between CORE and 
gCUT, because both techniques agree quantitatively for both restricted set of graphs. In 
our opinion there are basically two possible scenarios. First, an effective spin model can 
still be reliably derived at intermediate couplings using a hexagon expansion. Indeed, 
the minimal model only focusing on the leading one-hexagon graphs gives very nice 
results. It might then be possible that the full graph decomposition used for the gCUT 
by including all graphs up to seven sites does not give improved results, because the 
relevant low-energy physics is contained in graphs with longer loops (these are exactly 
the multi-hexagon graphs). Second, the trend seen in the full graph decomposition by 
gCUT is a real effect, i.e. a controlled derivation of an effective low-energy spin model 
for the QSL at intermediate couplings remains a big challenge. 

5. Conclusions 

Summarizing, we have derived effective spin Hamiltonians describing the low-energy 
physics of the Hubbard model on the honeycomb lattice [Eq. (IT])] at half-filling, that 
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Figure 10. (Color online) (a) Double occupancy per site obtained obtained with 
CORE effective model on finite clusters with N = 18 to TV = 30. (b) Derivative of 
the kinetic energy, as a function of t/U, obtained from the ground-state energy of the 
effective model [Eq. ^] on the same clusters. In both cases, comparison is made to 
ED (N = 18) and QMC (N = 72) jSTj data of the Hubbard model. 



has been recently shown to display a QSL phase [TTj . The effective models, obtained 
from the non-perturbative gCUT [20J and CORE [3(1 |3T] methods, are in good mutual 
agreement and seemingly well converged for the intermediate values of t/U yielding a 
QSL [11]. We find that the effective next-NN two-spin frustrating exchange J2 is not 
sizeable enough to destabilize Neel order [T51 (THl [LH Q2J [T9] and therefore cannot account 
for the existence of the QSL phase for Eq. 0. Instead, we find that the dominant sub- 
leading terms in the effective model are six-spin exchanges (see Fig. [2]), that may account 
for the emergence of the QSL behavior. 

Additionally, we have taken some first steps in trying to characterize the so-obtained 
effective model, by performing exact diagonalizations on small clusters. While the so- 
called Anderson tower of states is consistent with the occurrence of Neel order in the 
strongly interacting limit of t/U — > 0, the situation is less clear for t/U = 0.25 (that 
according to Ref. pj] yields a QSL), possibly an indication of a magnetically disordered 
state. 

Naturally, it would be desirable to unambiguously show that the herein derived 
effective model displays a QSL as its ground-state and possesses a gap to spin excitations. 
However, the original results for the Hubbard model [TTJ indicate that such spin gap 
is very small and thus one presumably needs to consider clusters comprising several 
hundred sites, so to surpass the spin correlation length and convincingly establish the 
nature of the ground-state. 

In face of this limitation, we can only conclude that the minimal effective spin model 
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derived in the present work gives very promising results, but that a comprehensive 
characterization of its ground-state is still missing. Accordingly, we hope that our 
results may stimulate further work in trying to understand effects due to the six- 
spin interactions Li, L 2 , ■ ■ ■ L 5 depicted in Fig. |2j One possible strategy would be to 
deform the here obtained effective model and to study its phase diagram in an enlarged 
parameter space, in the hope that the QSL would be further stabilized and the spin 
correlation length would be reduced within reachable system sizes. 
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